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Abstract 



Since the initial development of one-dimensional electron gases (1DEG) two 
decades ago, there has been intense interest in both the fundamental physics 
and the potential applications — including quantum computation — of these 
quantum transport systems. While experimental measurements of IDEGs 
reveal the conductance through a system, they do not probe critical other 
aspects of the underlying physics, including energy eigenstate distribution, 
magnetic field effects, and band structure. These are better accessed by 
theoretical modeling, especially modeling of the energy and wavefunction 
distribution across a system: the local density of states (DOS). 

In this thesis, a numerical Green's function model of the local DOS in a 
1DEG has been developed and implemented. The model uses an iterative 
method in a discrete lattice to calculate Green's functions by vertical slice 
across a 1DEG. The numerical model is adaptable to arbitrary surface gate 
geometry and arbitrary finite magnetic field conditions. When compared 
with exact analytical results for the local DOS, waveband structure, and 
real band structure, the model returned very accurate results. In zero mag- 
netic field, the local DOS plots from the model behaved as anticipated by 
theory; under a finite magnetic field, depopulation and waveband separa- 
tion were present in the model, also, precisely as was expected. The model 
was also used to investigate imaginary band structure and gave interest- 
ing results warranting further investigation. A second numerical model was 
also developed that measured the transmission and reflection coefficients 
through the quantum system based on the Landauer-Biittiker formalism. 
The combination of the local DOS model with the transmission coefficients 
model was applied to two current research topics: antidot behavior and 
zero-dimensional to one-dimensional tunneling. These models can be fur- 
ther applied to investigate a wide range of quantum transport phenomena. 
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Chapter 1 

Introduction 



The confinement of electron motion to a single dimension, an experimental 
achievement that opened new realms of physics, has garnered great inter- 
est since its development two decades ago [l j . As the limiting case in which 
current can be carried, one dimensional systems not only present fundamen- 
tal physics challenges but also portend numerous opportunities for practical 
application [2] . Most notably, they show potential to serve as the backbone 
of a future quantum-information-computation system [3]. 

A one-dimensional electron system is obtained by modifying a two- 
dimensional electron gas (2DEG). A 2DEG is created within a semicon- 
ductor heterostructure, a stack of a few different semiconductors that takes 
advantage of band structure to achieve only one allowable energy level in the 
z-direction. Most commonly the one-dimensional modification of a 2DEG 
is generated by the use of the easily-adaptable split-gate device [H [3] . Pio- 
neered in the 1980's, a split-gate is a strip of metal with a narrow slit across 



its width that sits atop the heterostructure (see Figure 1.1 ) |4J. By applying 
a voltage to the gate, the heterostructure beneath it is depleted, leaving 
only a tiny channel through which the electrons can move: that is, the area 
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split-gate electrodes 




2D electron 
system \ 



Figure 1.1: Split Gate Device. The metal gates sitting atop the semiconduc- 
tor hetero structure deplete the region underneath, leaving only a tiny channel 
in the middle in which electrons can be present: a ID channel 

underneath the slit. Since the width of this channel is roughly equal to the 
electron's wavelength, a one- dimensional electron gas (1DEG) is created. 

When van Wees et al [5] and Wharam et al [6] independently developed 
the first one-dimensional systems in 1988, their demonstration of quantized 
conductance resolved a three decades-old theoretical debate (see [7]) and 
sparked numerous novel investigations. Resulting studies in ID systems 
included investigations of the effects of high magnetic field and Zeeman 
splitting, of magnetic depopulation, of electric depopulation, and of tunnel- 
ing between ID channels [H El El llOj . These were followed by analyses of 
the behavior and number of occupied subbands and of the possible subband 
energies in ID systems jl H 112]. Inquiries into interference effects have been 
carried out extensively for the last 20 years. 

Yet, the physics underlying these 1DEG systems remains incompletely 
understood. |13|, I14j. Curiosity about IDEGs has only accelerated in recent 
years, fueled by discoveries such as spontaneous spin splitting (the 0.7 struc- 



CHAPTER 1. INTRODUCTION 



3 



ture) [15] . Further comprehension of experimental ID systems promises to 
be advanced by a deeper theoretical grasp of IDEGs underlying phenomena 

US]. 

What are these phenomena? They include the conductance, the capac- 
itance and the density of states [TBI HZ]- It is the density of states (DOS), 
in particular, that plays a fundamental role in understanding ID systems 
|18j . A measure of the number of states available in a given energy range 
per unit length (or area), the DOS can provide details about a system's 
wavefunctions, resonant states, thermodynamics, scattering amplitudes and 
transmission probabilities [HI [20] . 

Given its centrality to quantized transport systems, it is not surprising 
that much work has been carried out on the DOS in semiconductor het- 
erostructures. The total (bulk) DOS of a system has been investigated un- 
der periodic potentials [2T|, [22j [23] , under modulated magnetic fields [Ml EH], 
and in relation to localization length [26, 27 . The effects of spin [28] and 
non-resonant laser light [29] on the DOS have also been pursued. The two- 
dimensional DOS has been experimentally probed several times |30| I3"T1 132] . 
as has the tunneling DOS [33]. Countless other methods for calculating 
the bulk density of states for systems have been carried out, tailoring each 
model to meet specific material, dimensional or disorder constraints (see 
PS ESI QIJ EH E5l ES] ). 

The investigations listed above, however, have rarely focused on the local 
density of states. The bulk density of states measures the DOS averaged 
out across the sample; it returns a single value for every Fermi energy in- 
put. By contrast, the local density of states calculates the density of states 
independently at every given lattice point in the sample, returning thou- 
sands of values (each tied to a specific location) for every Fermi energy 
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input. Researchers have oft preferred studying the bulk behavior, because, 
in the words of one group of authors, they were simply "not interested in 
the details of the density of states." [21]. Yet, the local DOS is essential to 
understanding ID systems PJ. 

There have been some local density of states calculations of note (see 
|37|). Most relevantly, in 2003, Meyer et al carried out the first measure- 
ment of the local density of states in an extended ID system [2]. They 
measured the local DOS across a slice of the sample (a slice of a system will 
be more fully defined in subsequent chapters). Meyer then compared the- 
oretical predictions for the local DOS based on single-particle calculations 
with measured values of the local DOS. To their surprise, for their mea- 
sured values, they "did not find significant deviations from the calculation." 
Theoretical studies of the local DOS are therefore seen to provide a poten- 
tially rich source of both accurate and fundamental information regarding 
ID electron systems. 

This thesis will investigate the behavior of the one-dimensional, local 
density of states using a numerical approach. The work contained herein 
will concern itself not merely with the local DOS across the ID channel (as 
Meyer did), but simultaneously with the local DOS along the channel (in 
the direction of transport). Experiments on low-dimensional structures can 
only give conductance measurements, so one has to work backwards from 
these results in order to grasp the underlying physics. Yet, a numerical 
model — one that can compare conductance values and local density of states 
plots — provides a window into the structure of the system: into scattering 
effects, transport, magnetic field response, and imaginary band structure. 
In short, the local DOS gives a highly detailed portrait of transport (see 



Figure 1.2). With the aim of advancing knowledge about ID systems, this 
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DOS Across Slices 




x (lattice slices) 

Figure 1.2: Two-Dimensional Local Density of States plot for a 1DEG under 
the influence of a central surface gate. This plot is calculated using the 
numerical Green's Function model presented in this thesis. The system has 
width 101 lattice points and length 200 lattice points and is at Fermi Energy 
11 meV and field B =1 T. Black represents high density of states, white low. 
Depletion under the gate and the single subband that manages to pass along 
the edge are clearly visible. 

thesis derives and presents numerical models for finding both the local DOS 
and the transmission coefficients. The results confirm theoretical predictions 
with a high degree of accuracy. The model can be deployed to aid numerous 
inquiries and applications. 

This thesis contains four parts. Chapter Two provides an introduction to 
transport in 1DEG systems. In particular, the Landuaer-Biittiker formalism 
is outlined and so is the unifying work of Baranger and Stone. The third 
chapter describes the numerical method to be implemented, a method based 
chiefly on a technique forged by MacKinnon using Green's functions |38j . In 
the fourth chapter, a detailed description of the local DOS and transmission 
programs that are the products of this thesis will be presented. It is shown 
that the results from these numerical models match expected theoretical 
values for changes in magnetic field and gate voltages with robust accuracy. 
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Lastly, in Chapter Five, a selection of applications of the model is offered. 



Chapter 2 



Quantum Transport Theory 

2.1 Introduction 

To develop the numerical models presented later in this thesis, especially the 
transmission coefficients programs, a theoretical overview of quantum trans- 
port theory is required. This chapter begins with the basics of transmission 
in a perfect ID channel, then advances to the two-probe Landauer conduc- 
tance formula that can calculate transmission even in the presence of imper- 
fections in the system. From this two-probe form, the multi-probe Landauer- 
Biittiker formula is derived. The chapter then focusses on the work of 
Baranger and Stone who demonstrated the equivalence of the Landauer- 
Biittiker scattering formalism and the exact eigenstate (and Green's func- 
tions) formalism. Their work leads this thesis into a discussion of Green's 
functions (Chapter Three) which form the basis of the numerical model 
(Chapter Four). 
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2.2 Perfect ID Transport 

Classically, electron transport in metals was described by the Drude Con- 
ductivity. This equation, in which n s is the sheet carrier density, m* is the 
electron's effective mass, r e is the mean free time, and \i e is the electron's 
mobility, is given by: 

n s e 2 T e 

a = — = ensile (2.1) 

m* 

Though this equation works reasonably well in describing the transport 
for a 2DEG, it is not as effective for a 1DEG. The reason for this difference is 
that the dimensions of a 2DEG are typically much larger than the mean free 
path of the electron, meaning that one can use the average quantity r e with 
reasonable accuracy. Not so, however, for a 1DEG, where the width — and 
often the length, too — of the channel is typically less than the electron's 
mean free path. Since the electrons are usually confined electrostatically 
(e.g. by split gates), little to no scattering takes place along the width of 
the system. Longitudinal momentum is conserved. 

Since the Drude conductivity cannot accurately describe conductance 



in a 1DEG, a new formalism is required (see section 2.3). Before deriving 
this general new method, however, it is necessary to develop a description 
of basic transport in a perfect ID system by focussing on current flows in 
each direction. Assuming a perfect ID channel, the current flowing in one 
subband, i, does not scatter into any other, but flows cleanly through the 
system at its subband energy. If one considers a system with a chemical 
potential greater on its right-hand contact than its left-hand contact by an 
energy eV, then the current flowing to the right in channel % is given by: 
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(171 ' 

dl+ = -evlf(e + eV)^de (2.2) 

within an energy range de. v\ is the group velocity at energy e, /(e + eV) 
is the Fermi-Dirac function which gives the probability of an electron being 
in a given state, and ^ is the density of states, again at energy level e. 
One of the crucial features of this expression is that since v\ ~ and 
^Se ~ sir' S rou P velocity and density of states terms drop out. The 
system's dependence on energy, subband number, and momentum vanish 
from the expression for current, leaving only the Fermi-Dirac distribution 
function and constants: 

dlf = -^f(e + eV)de (2.3) 

Meanwhile, by an identical process, current flowing from the right-hand 
contact to the left-hand contact is simply: 

dl~ = ~lf(e)de (2.4) 

The net current flowing through energy subband i is the difference be- 
tween the right-flowing and left-flowing currents integrated across the full 
energy range: 

Ii=lJ™Jf(e)-f(e + eV))de (2.5) 
Since the Fermi-Dirac function is given by: 



f(e) = (2.6) 

e kT + 1 

where ji is the chemical potential, k is Boltzmann's constant, and T is tern- 
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perature, the function becomes a right-angled step function when the tem- 
perature is taken to zero, yielding the final transmission result: 

e 2 

h = jV (2.7) 
For N occupied subbands in the ID system, the conductance is: 

dl e 2 

g = w = n j < 2 - 8 > 

The above result is a fundamental theoretical insight into electronic be- 
havior. It does not, however, provide detailed information about transport 
beyond a perfect ID system. A lucid and simple formalism adaptable to 
studying wide-ranging quantum electronic systems including IDEGs is the 
Landauer-Biittiker method presented below. 



2.3 The Landauer Method 

The novel insight that earned Landauer his eponymous formula was the 
possibility of recasting a conductance problem as a scattering one |39l 140] . 
Instead of concentrating on the effect of applied electric fields, Landauer 
zeroed in on the transmission and reflection coefficients across channels in a 
quantum system [7]. If one considers a pair of ID electron systems, attached 
on either side to perfect Ohmic contacts and with an arbitrary potential 



region lying in between (see Figure 2.1) the two-probe Landauer formula 



can be derived straightforwardly. As was the case for the transport system 
analyzed above, a chemical potential on the right-side contact is greater by 
energy eV than the left-side contact. The current amplitude crossing the 
system to the right is positive and denoted by a + and b + , while the current 
amplitude moving to the left is negative and given by b~ and a~ . The total 
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Figure 2.1: Quantum channel with leads attached. Positive current moves 
right, negative current left, under the influence of the left lead's chemical 
potential, the right lead's chemical potential, (j,r, and passing through a 
quantum transport region with arbitrary effective potential, V(x,y). 

current amplitude that flows into the right contact, b + , and left contact, a~, 
can be written in matrix form: 



( b+ )- 




r+\ 




\ 


{ a- ) 


{ r 






J 



where t + , r + , and r~ are the transmission and reflection matrices. 

One can see that b + is simply the sum of transmitted portion of a + 
and the reflected portion of b~ , while a~ is simply the sum of the reflected 
portion of a + and the transmitted portion of b~ . Each of these ampli- 
tudes, a + ,a~,b + ,b~, are themselves vectors, their entries being, for exam- 
ple, af, 0,2,0,3, which are the amplitudes of a single subband i =1, 2, and 
3. Current transmitted to the right, in a narrow energy range de, is then 
determined by: 



dl+ = {a + H + h + a+) (2.10) 

with the above expression averaged over time. This expression represents 
the current that has crossed the system. The time average of (\af\ 2 ) is 
simply the familiar — | /(e + eV)de, and therefore the equation for dl + is 
given by: 
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dl + = -^Tr[t + U + ]f{e + eV)de (2.11) 

The trace of t + H + is calculated in order to sum the transmission co- 
efficients for each subband. Qualitatively, this sum makes sense, since the 
above equation is essentially calculating the portion of the current that is 
transmitted. By similar argument to above, for the right side of the system, 
one obtains: 

dr = -^Tr[t + U+]f(e)de (2.12) 

Using G = dl/dV, integrating for all energy levels, and once more taking 
the zero-temperature limit, one obtains the Landauer formula: 

G = ^Tr[t+tf+] = £ £ |tfcf (2.13) 

where Uj is the conductance coefficient for charge moving from subband i 
on the left-side portion of the quantum system into subband j on the right- 
side portion of the system. This derivation assumes an arbitrary effective 
potential and therefore applies to any ID system connected at both ends 
to ohmic contacts. In the ideal situation, these ohmic contacts behave like 
electron-emitting blackbody radiators. 

2.4 The Challenge of the Multi-Probe Formalism 

Landauer's initial 1957 paper determined the conductance formula to be: 



(2.14) 
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For decades this formula stood untouched [7j. It was only when Anderson 
jH], in 1980, and others 021 031 H SHJ S6] attempted to derive a multi- 
channel version of Landauer's equation that a conductance equation with 



G proportional only to T (2.13) and not G ~ T/R was derived |42[ 03]. 



Part of the concern with formulas like (2.13) had been that it gave a fi- 
nite resistance even for a perfect conductor. This finite resistance proved to 
be the result of the contact resistance: a resistance resulting from electrons 
with momentum distributions that did not match the quantized levels al- 
lowed in the channel trying to enter the ID system |47[ |4"8] . The attempt to 
develop a multi-channel formalism also grappled with the issue of whether 
the conductance measured at fixed current was a measure of the chemical 
potential at a set of charge reservoirs (source and sink far away from the 
system) or at an "effective chemical potential" inside the sample [7]. Ulti- 
mately, it was agreed that despite transmission being measured within the 
channel, it was acceptable to use the chemical potentials of the reservoirs 
outside the channel. This idea was demonstrated by Biittiker and is derived 
below E9|. 



2.5 The Landauer-Biittiker Formalism 

Biittiker put forward the successful multi-probe version of the Landauer 
theory by treating current and voltage terminals in a four-point probe set- 
up equally [H5J ED] • Swapping current and voltage probes, he demonstrated 
that such systems obeyed Onsager's relations |51| [52] (see jl9] for derivation 
of Onsager's relations). 

Biittiker considered four reservoirs, each at different chemical potential 
[ii with a fifth chemical potential, less than or equal to the lowest of all 
four Since states with energy below /j,q are filled,they cannot contribute 
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any net current to the leads. Therefore, the only relevant energy range is 
Afii = m — no above fio- The current injected by reservoir i into the system 
is, as was derived above: 

/* = %V = JA W (2.15) 
h h 

Part of the current is also reflected back into the reservoir. The magnitude of 
the proportion of current reflected is determined by a reflection coefficient, 
Ru, and results in: 

g 

Ii(Re fleeted) = ~ J^Rii^H ( 2 -16) 

Lastly, current is also flowing into the reservoir from other reservoirs. 
For each reservoir besides i (i + l,i + 2, etc.), there will be a transmission 
current with magnitude given by a transmission coefficient, e.g. Tyi for 
current injected into the system from reservoir 2 that ends up in reservoir 
1. Thus the impact on reservoir 1 from three other leads is given by: 



^{Transmitted) = ~T {T\2^^2 + 7l 3 A^ 3 + T14A//4) (2.17) 



In general, the net current flowing out of a given lead is therefore determined 
by the multi-probe Landauer-Biittiker equation: 

h = ~ R a)Vi ~ E r «Mij (2-18) 

The ^0 terms cancel out because the coefficients sum to zero. This can 



be seen clearly when (2.18) is written out as an TV" x N matrix, where N is 



the number of reservoirs. For a 3 x 3 system this would look like: 
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h 
\ h J 



1 — -Rn —T\2 —T\z 
— T21 1 — R22 —T23 

—T31 — T32 1 — i?33 



\ // .. \ ( „ \\ 

Mo 



Mi 

VMsy 



Mo 



(2.19) 



Every column and every row of the 3x3 matrix sums to zero, 1 — Ra — 
J2j Tij = 0, because physically speaking if the energy in every lead were 
the same, no current would flow. Thus, when any row is multiplied by a 
vector with identical entries — the //o vector — the result will be zero. This 
method of eliminating ^0 obviated the need for a "local chemical potential," 
allowing potentials to be calculated at the reservoirs alone. More generally, 
for N occupied subbands in lead i, the equation takes the form: 



(Ni - Ru)m - T i?Mj 



(2.20) 



This is the multi-probe formalism for an arbitrary quantum system with 
an arbitrary number of leads. Knowing only the reservoir chemical poten- 
tials, the transmission and reflection coefficients, and the number of conduct- 
ing subbands, one can calculate the conductance in any quantum system. 
This result is essential to the numerical model used to find the transmission 
and reflection coefficients implemented in Chapter Four. 



2.6 The Baranger and Stone Unification 

The final key piece of transport theory is the work of Baranger and Stone 
|53| . They showed the equivalence of the two chief approaches to linear 
response theory (/ ~ V), which are also the two approaches required for 
this thesis's numerical models: the Landauer-Biittiker scattering formalism 
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Figure 2.2: Quantum system with arbitrary number of leads. The shaded 
region is the arbitrary- geometry area of quantum interaction. Image from 

m 



(transmission program) and the exact eigenstate, Kubo-Greenwood, Green's 
function formalism (local DOS and transmission programs). Baranger and 
Stone made no assumptions except current-conservation, time-reversal sym- 
metry, and the non-interacting electron model. They began by finding the 
exact eigenstate form of the conductance coefficients, g mn that solve the 
linear equation I m = J2 n g m nVn, where leads m and n inject current into 



a quantum system with an arbitrary number of leads, Nl (see Figure 2.2). 



Then, they related g mn to the Landauer-Biittiker transmission coefficients 

± mn- 

They, following Biittiker, shunned the idea of "an effective chemical po- 
tential" lying somewhere inside the conductor. This "effective" or "local 
chemical potential," which plagued previous research had been developed 
in order to overcome the issue of particles scattering into voltage probes. 
Yet, conductors are necessarily out of equilibrium if current is flowing, and 
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hence theories about equilibration of different channels at the Fermi energy 
tended to be arbitrary. Baranger and Stone concerned themselves solely 
with current injected into the system. All that was required were the ap- 
propriate boundary conditions: a reservoir that was in equilibrium at fixed 
potential fj,, that was large enough so its potential would be unchanged by 
an additional particle, that inelastically scattered (phase-randomized) any 
particle entering the reservoir before returning it to the system, and that 
had a boundary with the sample that caused no additional resistance. 

Their derivation of g mn was carried out not by using conductivity, a spa- 
tially average quantity, but with the conductivity response function, a(x, x'). 
Since small spatial-fluctuations can greatly impact mesoscopic systems, it is 
important that (j{x, x') is a spatially-varying quantity describing the current 
density. It is also a function of states both at and below the Fermi surface; 
consequently, Baranger and Stone proved that the transport current is only 
a Fermi surface entity. Moreover, also unlike previous research S3] their 
derivation was the first to apply in a magnetic field of arbitrary strength. 



The conductance coefficients g mn are identified as (see Figure 2.3): 



g mn = - I dy m j dy'jtm ■ a(x, x) • x n (2.21) 

Cm Cn 

where x n is the unit vector parallel to lead n, y n points perpendicular to the 



lead, and C n is the cross section of the lead. As such, equation (2.21 ) gives 
a physically intuitive expression for g mn : it is the flux of the conductivity- 
response function from lead n into lead m that passes through the cross 
sections of each lead, f c dy m and j c dy' n . The full calculation of cr(x, x') 
is a laborious task whose expressions can be found as exact eigenstates or 
Green's functions in equations (40) and (75), respectively, in |53j. 

The conductance coefficients, g mn , are then manipulated so they are 
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Quantum Sample 



Figure 2.3: Transition from quantum sample into lead n. x n points in the 
direction of the lead n; y n is transverse to it. C n is the cross section of lead 
n, perpendicular to x n . Image from [53]. 

related to t mn , the subband transmission coefficients, by: 



where c is a subband in lead m, a is a subband in lead n, and —df/de is 
the derivative of the Fermi function. Taking the zero-temperature limit, one 
obtains: 



where T mn is the trace of the t + H + matrix. 

In this way, the exact eigenstate solution of linear-response theory stem- 
ming from a(x,x') and the scattering formalism derived from transmission 
and reflection coefficients are seen to be equivalent. Quantum transport cal- 
culations performed with exact eigenstates are identical to calculations per- 
formed with transmission and reflection coefficients. Each form has its own 






m 7^ n 



(2.23) 
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advantages. Green's Functions, which are derived from the exact-eigenstate 
formalism, are powerful numerical tools, while transmission coefficients give 
a physically intuitive description of transport. The method of using Green's 
functions for numerical simulations — the backbone of the programs of this 
thesis — is the subject of the next chapter. 



Chapter 3 



Numerical Green's Functions 



3.1 Introduction 

Green's functions are used to solve inhomogeneous differential equations and 
provide an effective method for analyzing the local density of states, con- 
ductance, and other transport-related properties of semiconductor systems. 
This chapter explores how they can be used to create a numerical model. It 
begins with an analysis of the discrete lattice (as opposed to continuous wave 
functions) and then considers the application of appropriate effective poten- 
tials to the system. This chapter then moves into an analysis of Green's 
functions: their definition and their utility, how they are used to solve a 
Hamiltonian system, how they are used iteratively to calculate transport 
properties, and lastly the appropriate boundary conditions. The iterative 
process, it should be emphasized, is crucial, as it allows one to calculate the 
energetics of the entire system (e.g. 100,000 lattice points) by repeatedly 
multiplying only vertical slices of the system (each with only 100 lattice 
points); this permits the multiplication of matrices of the order of 100 x 100 
rather than having to perform an inversion of a 10 5 x 10 5 matrix — a massive 
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numerical task. 

3.2 Creating A Discrete System 

3.2.1 Discretizing the Schrodinger Equation 

In order to create a matrix representation of a quantum system, the relevant 
Schrodinger equation must be discretized. To allow for maximum generality, 
a magnetic field is applied, perpendicular to a two-dimensional system, with 
magnetic vector potential, A. The Landau gauge, A = (—By, 0,0), is used. 
Employing the Peierls substitution of inserting eA into the Hamiltonian |54[ 
[55] , the Schrodinger equation for the two-dimensional system is determined 
by: 

2b(-' h l I +eB ^*-&p- +v * = E ' p (»» 

A lattice constant, a, is introduced. This transforms the equation into: 



' Q d + ieBya\ 2 
v dx h J 



d 



ip + 



2m*a 



*„2 



*„2 



dyj r ' h 2 



2m* a 

~h 2 



■E$ (3.2) 



The squared terms in (3.2) are replaced by a second order Taylor poly- 
nomial for e x + e~ x using x 2 ~ e x + e~ x — 2: 



/ „ 9 i~,y n 8 i-iy\ , / a d_ - a e_\ , 2m* a 2 T r , 2m* a 2 ,_, , 
4ip-[e a^e^ + e a ^e e s « + e a « J VH -3— V>/> = 2 £-0 

(3.3) 

where 7 = 

Now the system is made discrete by deploying a Taylor expansion once 
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again. Using the fact that, to second order: 



and that 



e^(z))^(*) + a^ + ^g (3.4) 



^ + a)«^(x)+a^ + |^| (3.5) 



one sees that e a sx (^(x)) ~ ip(x+a). The corresponding result e a sx (^(x)) 
ip(x — a) also holds. Replacing x and y with lattice points n and m related 



by x = na and y = ma (see Figure 3.1), the discretized Hamiltonian is 
obtained: 

e nm V>n+l,m + e~ n ">n-l,m + 1pn,m+l + 1pn,m-l + Vlpn,m = £lpn,m (3.6) 

where v = — 2 T° 2 V", e = 4 — and ip n +i 7 m represents the wavefunction 

one lattice point right of the wavefunction at point (n, m). Using standard 
error approximation methods for a Taylor series for e a e* , the error will be no 
greater than R n = ^-e^ where £ is less than a, and for e~ a 8* the error will 

3 

be no greater than i? n = Typically, a is of length 5 nm, making these 
errors very small indeed. As will be shown in Chapter Four, the numerical 
model's results and the results expected by theory match extremely closely. 

This discrete result is known as the tight-binding approximation. If one 
inserted Vn,m = e lkxCLn e lkyarn as a solution and set v = 0, one would end 
up with cosine bands as the dispersion curves: E = 2cos(fc r a) + 2cos(k y a). 
Thus the solutions are tightly bound in the curvature of the cosine bands. 



Having found a discrete form of the Hamiltonian (3.6), the system can now 
be translated into matrix form and the method of numerical Green's func- 
tions described. First, however, the value of the effective potential, v, is 
calculated. 
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Figure 3.1: Lattice used for computation. Each dot represents a lattice point. 
The numerical model is carried out by multiplying matrices representing ver- 
tical lattice slices (e.g. slice n). Positive current is taken to be flowing in 
the positive x- direction. 
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3.2.2 Calculation of the Effective Potential 

To calculate the effective potential in a quantum system, the most thor- 
ough mechanism would be to calculate the potential self-consistently. This, 
however, is not the method used here. Self-consistent potentials require the 
simultaneous solution of the Schrodinger and Possion charge distribution 
equations, a numerically intensive process. As such, though self-consistent 
potentials produce useful results (see, for example [56]), they do not well 
serve a model designed for rapid adaptation to a large range of geometries 
and surface gates in quantum systems. 

Instead, a powerful and easily malleable tool for calculating the potential 
can be found in the work of Davies et al |57j . Their model is specifically 
tailored to measuring the effects of gates placed on the surface of a 2DEG. 
They do not factor in the contribution of the fields generated by the electrons 
themselves, but their results are nonetheless very practical and accurate for 
numerous reasons detailed in |57] . 

Their model derives from the solution to Laplace's equation, V 2 = 0. 
The first boundary condition is that c/>(r,0), the potential, is equal to V g , 
the applied gate voltage, and the depth z = is the surface where gates are 
put down. The second condition is that d<j)/dz = in the limit z — ► oo. 
Then the two-dimensional Fourier transform is applied to cp(r, 0) turning 
it into </>(q, 0). Given that z must decay exponentially in order to satisfy 
dcp/dz = as z — > oo, the general expression for the transform is given 
by: 0(q, z) = 0(q, 0)e - ' 9 ^. This multiplication of the Fourier Transform, 
of course, is the same as convolution in real space. As a result, taking the 
inverse Fourier Transform yields: 
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This general equation can then be manipulated for a whole host of re- 
sults. The geometries of the gates simply need to be given in polar coordi- 
nates r = (r,9). In general, the results are arctan(x, ?/) functions, resulting 
from the integral across the surface of the 2DEG. 

The most common gate deployed in this thesis's calculations was a finite 
rectangular gate. Its effective potential is given by: 




= g{x-L,y-B)+g(x-L,T-y) (3.8) 



+ g(R-x,y-B) + g(R-x,T-y) (3.9) 

where g(i,j) = ^arctan(^) and R = \/i 2 + j 2 + d 2 . The depth of the 
2DEG below the surface is d, and L, R, B, and T are the values of the left, 
right, bottom, and top edges of the rectangular gate. Naturally, one could 
create an arbitrary number of gates and simply sum their effective potentials 
by the power of the superposition principle. 

While the Davies et al formulation allows for the calculation of proper- 
ties for gate designs of all varieties and geometries, certain calculations are 
best carried out with a confining potential free of surface gates. In such 
a scenario, the confining potential of an infinite square well or a simple 
harmonic oscillator can be used with effective results (the potentials being 
set up transverse to the current). Their numerical implementation is dis- 
cussed briefly in the next chapter. The saddle point potential — c/>(x, y) = 
(j)o — ^mu! 2 x 2 + ^mu> 2 y 2 — is another very effective model for the potential 
arising from a split-gate [58 . 
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The groundwork for the numerical method has been laid by developing a 
discretized quantum lattice and an effective potential. The Green's function 
numerical technique is now presented. 

3.3 The Green's Function Numerical Method 

3.3.1 A Green's Function Primer 
3.3.1.1 The Definition 

Green's functions are implemented to solve inhomogeneous differential equa- 
tions. Consider a partial differential equation of the form: 



where £ is a linear operator on y{r) and T is the inhomogeneity. The 
solution, y(r), is written in terms of the Green's function, G(r,r'), and its 
product with the inhomogeneity: 



Hence, T{r) = / CG(r, r')T{r')dr' ', and as a consequence of the definition 
of the dirac-delta function: 



The Green's function can be thought of as the inverse of the linear op- 
erator, as is especially clear when the differential equation is cast in matrix 
form, CG = I, where I is the identity matrix. 



Cy(r) = F{r) 



(3.10) 




(3.11) 



CG(r, r') = 5{r - r') 



(3.12) 



CHAPTER 3. NUMERICAL GREEN'S FUNCTIONS 



27 



3.3.1.2 The Eigenvalue Equation and s 



Equation (3.6) can be recast as the energy eigenvalue equation Hip = eip. 



This equation can be rewritten (e — H)ip = 0, and in this form, e — H plays 
the role of C. The Green's function is therefore given by: 



(e-H)G(r,r') = 5{r - r') 



(3.13) 



Normally, however, a complex energy z = e + is is defined and used 
to replace e. s is made infinitesimally small, in order to avoid having an 
impact on the numerical result |38j . The reason for the inclusion of s can 



be understood when 3.13 is rearranged for G(r,r'): 



G{r,r',z) = ——jj.6(r-r') 

~ \ Z _ R (3-15) 
= Y, Mrmr ' ] (3.16) 



Where the substitution of e n for H in (3.16 ) is made because H\(f) n ) = e n \4> n ). 



As can be seen, since H is Hermitian and therefore e n is real, G(r, r') is 
analytic everywhere except at the eigenvalues of H. G(r, r') has poles at 
these discrete eigenvalues, e n . Consequently, to avoid the problem of having 
to calculate residues wherever z = e n , s is added to e n . In calculations 
performed in Chapter Four, s was set to 1CP 18 : small enough to have no 
impact on the results, but large enough to avoid computational error. 
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3.3.1.3 The Dyson Equation 

One of the most important properties of Green's functions is their simple re- 
formulation when a perturbation is added. Consider a system with solution 
Go = z ^h = ( z ~ Hq)^ 1 . A perturbation with Hamiltonian H\ is added. 
The total Hamiltonian is now given by H = Hq + H\ and the total Green's 
function is: 



G = {z-H Q -H 1 )- 1 = {G 1 -H 1 )- 1 (3.17) 

Multiplying both sides by the inverse of the right hand side and then by 
Go one arrives at G — GqH\G = Go, or alternatively: 

G = Go + GoHiG (3.18) 

This relation is the well-known Dyson equation. Since adding a new lattice 
slice to the quantum system is adding a new perturbation, H\, the Dyson 
equation plays a fundamental role in the iterative method described in sec- 
tion 13.3.31 

3.3.2 The Green's Function and the Hamiltonian 

The Hamiltonian of the system must account for every point in the N x M 



lattice (Figure 3.1 ). The total Green's function matrix for the whole system 
must therefore be of the same dimensions. Fortunately, calculating these 
enormous matrices is not required. Rather, if the system is divided up 
into N vertical slices, each slice having M lattice points, then the relevant 
matrices are the matrices of each slice, of dimension M x M instead of 
NM x NM. Inverting an NM x NM matrix to solve for every entry of G 
would be an enormous numerical calculation for even a modest system of 
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./V = 500 and M = 100. Instead, the matrices for the Hamiltonian H, the 
energy Z, and the Green's functions, G are all made to correspond to only 
a single vertical slice. The matrix relationship relating all the lattice points 



on slice i to all the lattice points on slice j (see Fig. 3.1 ) is given by: 



[Z — H^j] Gjj 



H 



Gj+ij — Hj^-iG 



(3.19) 



Given the nearest neighbor approximation form of 3.6 in which only 



the effects of neighboring lattice points impact the Schrodinger equation for 



that lattice point, 3.19 is appropriate here. As a result, only three terms 
are present above: z — H on the slice, H one slice to the right, and H one 
slice to the left. For clarity, it helps to write out [Z — H^j]. On slice re, with 
effective potential at point m given by v n ^ m , [Z — H^,] is: 



Z~V nt _M +1 "I 



-1 Z - U n n -1 



4 Z ~"n,f 



(3.20) 



The above can be thought of as the matrix representation form of (3.6), 



where each row of the matrix corresponds to the Schrodinger equation for 
a particular lattice point located at position (n, m) in the lattice. The 
e^^Vn+i.m and e _l7m V ; n-i,m terms are reserved for the Hj^+i and Hj^-i 



matrices, respectively. There are a few features worth noting in (3.20). First, 



e has been replaced by complex z. Second, the off-diagonal -1 terms rep- 
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resent, reading across a row, the if) n ,m—l an d VVim+i terms of Hamiltonian. 
Thirdly, the effective potential is centered around zero, from lattice position 



-4f to 4f, instead of from to M. 



The reason for centering each slice's label around zero is to make the 
matrices symmetric as is evident in the forms of the other two Hamiltonian 



matrices in (3.19). The Hamiltonian linking one slice to its neighbor on the 
right is given by: 



Hj i-f-i — V 



-i 7 f 








e H(-f +1) 



(3.21) 



where, as before, 7 



eBa 



-. The third contribution to (3.19) is defined by 



H 



Vt. In both V and V* the matrix is ordered from 



M 



to 



M 



for 



symmetry. Of course, adding a small translational shift, to the magnetic 
vector potential A = (—B(y + (),0,0), would not change the magnetic field: 
B = Bz. 

A final critical feature of these matrices is their translational invariance. 
That is, except for the varying effective potential, every one of the three 
matrices [Z — H,^], V, and V^, is identical no matter what slice i is being 



calculated. This, of course, is logical given that (3.6) has an identical form 



for every lattice point and (3.19) has an identical form for every slice. In the 



presence of a non-translationally invariant potential, however, neither H^j 
nor the Green's function matrices themselves will be the same for every slice. 
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This is an important result, because otherwise the Green's functions would 
provide no information about the energetic changes across the system. 

3.3.3 The Iterative Process 

To calculate the Green's functions across a sample an iterative process is 
used. First, an initial matrix, Go o> is determined using boundary conditions 



(see section 3.3.4), and then each ensuing Green's function, Gi.i, G2,2 etc., 
is calculated from the previous Green's function. As explained earlier, this 
iterative process derives from the Dyson equation. For the purposes of the 
MacKinnon method, this equation should be recast: 

Gg +1) = + Gg for (i, j < N) (3.22) 

Where as before, represents the interaction between slices i and j, and 
superscripts n and n + 1 represent the number of slices incorporated into 

(3) 

the calculation thus far. Thus, for example, G 2 3 is the Green's function 
representing the interaction between slices 2 and 3, calculated after iterating 
to slice 3. It is worth noting that this equation is very well-behaved upon 
repeated application. Numerous tests carried out during the writing of this 
thesis's local DOS program consistently showed that the Green's functions 
converged as iterations were carried out for systems of various lengths. 

Generally, each of the applications of Green's functions require iterating 
from slice to slice N. The goal is therefore to be able to obtain any 
Green's function in the system having iterated to slice N, since this takes 
into account all the energetics of the system. For example, > is the 
Green's function linking slice 2 to slice 3 after the iteration has carried all 
the way through to the edge of the system at slice N. It is the definitive 
value for that interaction (leaving aside boundary conditions), as opposed 
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to G5> 3, which is a matrix that does not take into account the impact on 
the slice 2-3 interaction resulting from the Hamiltonians of slices 4 all the 
way through N. Gg 3 differs from G^ 3 ] for non-zero field, because on each 
iteration, every MxM Green's function matrix to be calculated is multiplied 
by either V or , and these perturbation matrices equal the identity matrix 
only at B=0. 

The details of carrying out this iterative process to slice N for the density 
of states and transmission coefficients calculations require substantial spe- 
cific explanations. Consequently, they are put off for Chapter Four. Here the 
general iterative method, moving on from the Dyson equation, is explicated. 

The Dyson equation, though very useful, is recast into four equations 
derived from it [3"H] : 



,(n+l) 

'n+l, n+l 

n (n+l) 
"i,n+l 



G 



(n+l) 
n+l,j 



[Z-H^-VtGgV]- 1 

r> ( n ) 1 r 1 ( n \m ( n+1 ) \r\n ( n ) 
4 + „ v^ n+lin+1 v 



1,3 



i,n 



GSvGiti^i (i < N) 
GSS^VtGg (j < N) 



(3.23) 

(i,j<N) (3.24) 
(3.25) 
(3.26) 



These four equations divide up the total Green's matrix that takes into 
account n + l slices, G^ n+1 \ into four regions. The total Green's functions 
matrix at this iteration isann + lxn + 1 matrix with each entry itself being 



an M x M matrix. Each equation (3.23)-(3.26) is capable of determining 



only certain MxM matrices in the total matrix, though together they can 



find them all (see Figure 3.2). 



Equation (3.23), finds the self-interaction energy of the n + l slice, the 



last row, last column entry of total matrix after n + l iterations, G^ n+1 ^. 
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Figure 3.2: The Range of the Four Iterative Green's Functions. The figure 
on the left describes which equations (of (3.23 3.26)) are used to calculate 
the entries of G^ from G^ . The right figure describes which equations 
are used to calculate the entries of G^ from G\ N ~ l > . Since each Green's 
function matrix, g\^\ is an M x M matrix, this right figure represents the 
total Green's function matrix, G^ N \ an NM x NM matrix. Only certain 
of the Green's function matrices within this total matrix are needed for most 
numerical applications; thus, the iterative method ends up saving a great deal 
of calculational time. 
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Equation (3.24) is capable of finding any entry in the total matrix except for 
the last row and the last column. It cannot find them all at once, however, 
and specific values for both i and j must be implemented in order for the 
iterative process to work. As will be shown in Chapter Four, this is not 
problematic, as the Green's functions for the density of states are almost 



always sought for the case of slice i = j only. Equation (3.25) calculates the 



final n + 1 column of the total matrix, while (3.26) finds the n + 1 row of 



the total matrix. The derivation of each of these equations from the Dyson 



equation is offered in Appendix A. From these four matrix relations (3.23) 



(3.26) every Green's function relating any two slices of the quantum system 



can be found by iteration. 

3.3.4 Adding Leads: Green's Function Boundary Conditions 

The issue of the boundary conditions remains. Though the iterative process 
allows one to calculate every value of the Green's function across a system, 
it does not take into account the interactions of the system's edges with the 
leads. The most relevant boundary conditions to consider are where two 
semi-infinite metal leads are attached to each end of the quantum system 



(see Figure 2.1). These conditions are imposed below and were those used 



in the numerical models of Chapter Four. 

3.3.4.1 Determining and Sorting Eigenvalues and Eigenvectors 

To begin applying boundary conditions, the matrix form of the Green's 



function equation (3.19) is recast as an eigenvalue problem: 
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\ 






V 




V 






/ 



VtGi_i 



V 





(3.27) 

To solve for G, the homogenous case of the same eigenvalue problem is 
considered. For eigenvalues a and eigenvector matrices U a and U b one 
obtains: 










(3.28) 



A single matrix whose eigensolutions are sought must be formed. To 
accomplish this, both sides are multiplied by the inverse of V. Using the 
fact that V -1 = V"!", the result, known as the transfer matrix, is obtained: 



( Vt(Z - H) -Vt ^ 

vt 



(3.29) 



The eigenvectors contained within U a represent the wavefunctions of 
the quantum system. They can be sorted according to the magnitude of 
their corresponding eigenvalues. If a > 1, then the wavefunction is an 
evanescent mode traveling with positive momentum, while those vectors for 
which a < 1 are evanescent modes traveling with negative momentum. In 
the case of a = 1, the corresponding eigenvector is a conducting mode. 
These current carrying modes must be normalized. 

Eigenvector U a is separated into two M x M matrices, U + and U_, 
corresponding to the direction of the wavefunctions' momenta. Eigenvalues 
a are divided up appropriately into a + (for a > 1 and positive momentum 
values of a = 1) and a_ (a < 1 and negative momentum values of a = 
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1). The current-carrying modes, i.e. those wavefunctions with a = 1, are 
determined to have positive or negative momentum using the definition of 
current in a quantum system. If 2ip^ m Im(aV m ) > 0, then the current is 
positive (moves right). If this expression is < 0, the current is negative 
(moves left). 

3.3.4.2 Applying Boundaries 

With a+, U+, a_, and U_ in hand, an expression for the Green's functions 
corresponding to each boundary can be derived. Ggj]° is defined as the 
Green's function from the end of the semi-infinite lead on the left side of 
the sample to the zeroth slice of the sample. G~^°^ represents the Green's 
function from the final slice N of the sample to the end of the semi-infinite 
lead on the right side of the sample. 

The left lead, Gq^°, is considered first. The boundary condition here is 
that the Green's function relating slice — oo to slice must go to zero. The 
reason for this is that the Green's functions must decay into the lead: the 
system's energy ought to go to zero as one moves infinitely far away from 
the quantum sample and into the current injector. From the comparison 



of the inhomogeneous and homogeneous eigenvalue formulations, (3.27) and 



(3.28 ), and the fact that according to the Bloch theorem for a regular lattice, 



Gj+ij = aGj j, it is identified that: 



Gr°° = U+a^A (3.30) 

for two arbitrary horizontal positions in the lead, % and j, where i < j, and 
where A is a matrix of coefficients. \a\ > 1 in order to satisfy the boundary 
condition (given that i < j). To derive an eigenvalue expression for G7j° 
independent of A, two cases of different initial values i = = and 
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0, j = are considered and plugged into (3.27). 



VG ,o = (Z-H)G_ li0 -VG„ 2 , -0 (3.31) 
VGi,o = (Z - H)G 0)0 - VG_ lj0 - 1 (3.32) 



Immediately, VGi o goes to zero because i > j. The values for Goo, 



G_io, and G_2,o are determined by use of (3.30): 



Go,o 
G-i,o 

G_2,0 



U+A 



U+q7A 



U+a7 2 A 



(3.33) 
(3.34) 
(3.35) 



Substituting these terms in and (Z — H) out, simplification leads to: 



A = a7 1 U+- 1 V t (3.36) 



And thus: 



(3.37) 



A nearly identical process is applied to calculate the Green's function in 
the right lead. The Green's function again must decay into the lead, but for 
this to hold true here i > j and |q| < 1. Comparing the inhomogeneous and 
homogeneous equations again one obtains: Gf°° = U_crL l A. Substituting 
in two sets of i and j, the result is: 
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G+°^ = U_a_UZ x V 



(3.38) 



The expressions Gg q° and G~jy°jy are inserted into the appropriate places 
in the numerical process. Carrying out the iterative process, one first inserts 



'0,0 



into the right-hand side of (3.23) in the place where Gqq sits in (3.23). 



This will yield the result G^. The G o° expression is also inserted into 



(3.25) in the places where G-q^ sits, and is inserted into (3.26) in the place 



where Gg°j sits. Applying this expression in those cases gives the Green's 
functions, Gq 1 ] 1 and G^q, respectively. Furthermore, Gq™ must be inserted 



into (3.24) on the first iteration in place of G-^, G-° ' 1 , and Gq }, returning, 



Gg 1 ^. Applying Gq o° into these four equations on the first iteration, allows 
one to take into account the Green's functions running all the way into the 
left lead. 

The implementation of the right lead has one slight nuance. It is inserted 



only into (3.23). Since ( 3.24 )-( 3.26) each depend upon the G|^„ +1 matrix 



emerging from (3.23), implementing the right lead Green's function once in 



(3.23) is sufficient. The definition of the Green's function as (Z — H ra 



+1J 



is used and (G^°^) is inserted into its place: 



G&h W i = [(G^r 1 - VtG^V]- 1 (3.39) 

This completes the formalism of the Green's function numerical model. 
Iterating Green's functions provides a powerful and efficient tool for cal- 
culating fundamental quantum mechanical properties of electronic systems. 
This is evinced in Chapter Four, where numerical Green's functions are ap- 
plied to calculating the local density of states and transmission properties 
of one- dimensional quantum samples. 



Chapter 4 



The Density of States Model 

4.1 Introduction 

In this chapter, two programs for calculating properties of quantum systems 
are presented and their results are analyzed. The first program calculates 
the local Density of States (DOS) of a ID quantum sample, the second the 
transmission coefficients of a ID quantum sample. The local DOS is a fun- 
damental property of a quantum system and the use of numerical Green's 
functions in a discretized lattice provide an effective probe of its behavior. 
The two programs were written in Visual C# (C-Sharp), using Center- 
Space. Matrix auxiliary code and NPlot Graphics code. All code developed 
is the work of the author with two exceptions: a piece of code that made 
assigning eigenvectors to positive and negative momenta matrices more ef- 
ficient, and a piece of code in the transmission program that rapidly carried 
out a series of multiplications for the transmission and reflection coefficients. 
These two pieces of code were based on the work of C. H. W. Barnes [59 . 

This chapter will begin by discussing the local DOS and the method by 
which the program was compiled. It will then compare results generated 
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from the numerical model with expected theoretical results and will con- 
firm the model's very high degree of accuracy. Finally, the method of the 
transmission program will be presented and its results analyzed. 

4.2 The Density of States Program 
4.2.1 The Density of States Function 

The density of states is a function of energy that measures the number of 
states available in a given energy range per unit length. The number of 
available states depends upon the number of occupiable states in k-space, 
and the DOS is in general a measure of how closely packed energy levels 
(and their corresponding wavefunctions) are in a quantum system. The 
one-dimensional total density of states is given by: 



where, as before, N is the length of the ID system in lattice units, M is 



its width in lattice units, and i is a slice of the system (see Fig. 3.1). This 



equation (4.1) is the density of states averaged out across the sample. To 
determine the local density of states, the density of states is individually 
measured at every lattice point m on slice i, and the relation to the Green's 
functions becomes: 

PH = ^Im (G<f) (4.2) 

evaluated at entry (m, m) in the G matrix. 

This result emerges from an elegant physical and mathematical argu- 
ment highlighting the fundamental nature of the DOS. The density per 
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energy value e is J2 n ^( e ~ e n)i where e n is an eigenvalue of the Hamilto- 
nian. The DOS then is the product of this density per energy level and the 
corresponding wavefunction probability map, ^(r) 2 !: 

pi r) = J2S(e -e n )|V>(r) 2 | (4.3) 

n 

To link this to Green's functions, one must return to the original Green's 



function formalism. Due to the poles in the Green's function in (3.16), 
a branch cut of G(r,r',z) is taken along the real-axis, creating two new 
Green's functions, the retarded and advanced Green's functions: G + = 
lim s ^ + G(e + is) and G~ = lim s ^ + G(e — is). The retarded Green's func- 
tion, z = e + is, was selected for use in the calculation. Using the identity: 



lim = P(-) T m5(x) (4.4) 

»/^o+ x ±iy x 

where P is the Principle Value Term, one can recast the retarded Green's 
function. Defining x = e — e n and y = s (since — — = — ^-ttt) one arrives 
at: 



G + (r, r\e) = PY J - ^v^ ^ 5(e - e n )MrW n (r') (4.5) 

n 6 £ ™ n 

And therefore: 

p(r) = - Im (G+(r, r', e)) (4.6) 

This demonstrates the very close link between Green's functions and the 
DOS. By calculating the Green's function at every ijj nm in the lattice, the 
density of states is found. No simplifications or approximations whatsoever 
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(except to create the discrete lattice) need be made. 
4.2.2 The Method of the Program 

In creating the program, most of the input constants, variables, and matrices 
are simply defined and plugged in to meet the specifications of the system. 
The one area that is tricky, and which this section will spend substantial 
time dealing with is the iteration to the proper slice in order to calculate 
the local DOS. 

To begin, the lattice width M, lattice length N, lattice spacing constant 
a, and applied perpendicular magnetic field B are defined. Each of these 
variables can be readily varied. The effective mass used was that of GaAs, 
0.067e. The constants v = -^^V and e = 4 - ^£e as described before 
are used, where E is the input voltage in meV. Infinitely small complex value 
is is added to e. Matrices V, V^, and Z are initialized with appropriate 
values, as explained in Chapter Three. 

The Hamiltonian matrix is an effective potential matrix with -1 on the 
off diagonals. The effective potential values v run down the main diagonal. 
Recall that this matrix represents the Hamiltonian for a slice of the system 
only. For the case a translationally-invariant potential (e.g. infinite square 
well or harmonic oscillator), every slice will have the same Hamiltonian, each 
main diagonal entry corresponding to the potential at a lattice point as a 
function of m. For example, the first matrix entry corresponds to lattice 
point m = —M/2 at the bottom edge of the system, and therefore must 
reflect the potential — including the effect of the magnetic field — a distance 
M/2 from the system's center. 

Specifically, for the case of an infinite square well, v = all the way 
down the diagonal, while for a harmonic oscillator potential, V = ^(m — 
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M-l\2 



, giving a minimum at matrix entry, M/2, which corresponds to the 



center of the system, row m = (see Fig. 3.1). For the more complicated 
translationally- varying potential, like that emerging from a surface gate, the 
matrix must also be made a function of N, and it must be recalculated at 
every iteration. 

With the essential matrices in hand, the eigenvalues and eigenvectors of 
the transfer matrix are found and sorted by modulus. This is a numerically- 
direct but programatically-heavy technique. In short, code was written to 
effectively sort the N eigenvectors and eigenvalues by their modulus and 
keep the eigenvectors and eigenvalues appropriately paired. A separate piece 
of code deals with those eigenvalues whose modulus is equal to one — the 
current-carrying modes — and sorts them into positive and negative current. 
Still another piece of code sorts the a > 1 and positive momentum values 
of a = 1 into a + matrix and the a < 1 and negative momentum values of 
a = 1 into a_. After a+, a_, U+, and U_ have been found and sorted, 
and after they are used to determine Gqo°, the iterative process begins. 



Looking back to (4.2), one sees that in order to calculate the density 
of states iteratively, one must find the G^P matrix. This is the Green's 



function self-interaction at slice i, analogous to G + (r, r',e) of (4.5) with 
r = r'. Recall that the superscript in G^P indicates this Green's function is 
not the result of iterating to slice i (which would be G'fJ ) , but is an iteration 
all the way from slice to N. In fact, since boundary conditions are included, 
it is an iteration from — oo to oo. Therefore, to calculate the local density 
of states using the iterative method, a process is needed to calculate the 
Green's functions at any slice inside the system while accounting for an 
iteration all the way to oo. 



This method will naturally rely upon (3.24), with i = j and n = N; 
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however, finding this value depends upon iterating the other equations in the 
proper order. The method developed is as follows. The program calculates 



(3.23 ) from slice (using G Q ^° as a starting point as explained in 3.3.4 ) until 



a particular slice i = w. It is at slice w where the local density of states will 



be calculated. Note, again, that when calculating (3.23), if the potential 



varies with N, the Hamiltonian must be recalculated at every iteration, 
thereby slightly increasing calculation time. 

Once the program has iterated to slice w, the value of G-K, is stored 



and then (3.23) is applied again. Then (3.24) is applied, using the stored 



input of Gw)i) three times and G[ ) ^ t 1 1 ^ +1 to obtain G^u^. Subsequently, 



(3.25) and (3.26) are applied, employing G™ w and G^+i^^ once each to 

(w+l) 
w+l,w 



calculate G^w+i an< ^ ^w'+iw respectively. In the next iteration, applying 



(3.23) will yield G^^w+2- When iterating (3.24) again, it will be given by: 



C(w+2) _ r>(w+l) , pi(uH-l) vf( (t»+2) vtr , ( M,+1 ) (A 7 s ) 

^w,w — ^*w,w "T ^w,w+l v y *w+2,w+2 v ^w+l,w \*-') 

This process explains why finding Gw,w^ depends upon the other iter- 



ative relations. It requires (3.23) and its G^+2w+2 result, (3.25) and its 



G-ww+l resm t, (3.26) and its G^J 1 ^ result, and (3.24) itself and its G^™^ 



result. Equation (3.24) must be iterated after (3.23), but before (3.25), 



and (3.26), since (3.24) relies upon their values from the previous iteration. 



These four equations are then applied over and over again until slice N, at 
which point the desired G^J, matrix is obtained. 

At this point the Green's function in the right lead, G^y°^-, is inserted into 



(3.23), according to the procedure described previously; then this result is 



used to find (3.24). This yields the ultimate result: G+ 1 ??, which can also be 



thought of as Gww \ since it takes into account the Green's functions 



in the leads and in every slice of the sample in between, running from — oo 
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to +00. Each diagonal entry in is the value of the Green's function at 

lattice point to, beginning from the bottom of the lattice slice and moving to 
the top. Taking the imaginary portion of this Green's function according to 



(4.6) yields the local DOS at that lattice point alone. In order to calculate 
the density of states in every slice of the sample, the code found the local 
DOS for slice w = 1, then started over again and calculated it for slice 
w = 2, repeating the process until slice N. Thus the local density of states 
across the entire ID system, in both x and y, was determined. The code 
stored every lattice point's local DOS in a two-dimensional array where each 
column corresponded to a slice; these values were then fed into a graphics 
compiler. 

Given this method, the degree to which N is increased or decreased 
will lengthen or shorten the calculation time by the fractional change in N 
squared. For example, the system doubled in length from 200 to 400 lattice 
slices, each iteration of a slice would require iterating twice as far to get to N 
and then the program would also require twice as many iterations in order 
to include all ./V slices. There is, unfortunately, no way around this fact, as 
the density of states not only must be calculated at each slice but also must 
take into account — on every iteration — the Green's function linking every 
slice to its neighbors. A typical calculation, with M = 41 and N = 200, 
took five minutes on a standard desktop computer. 



4.3 Results from the Density of States Program 
4.3.1 Testing the program's accuracy 

The Green's function numerical model developed for this thesis matched 
expected results for the local DOS with extreme precision. First, the shape 
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and response of the wavefunctions to varying magnetic field were examined. 
Next, the properties and behavior of the real bands of the system at varying 
magnetic field were considered. Lastly, Green's functions plots for a slice 
of the local density of states were compared with results using the analytic 
formula for the DOS derived from the carrier density. In all three cases, the 
Green's function model demonstrated great precision and reliability. 

4.3.1.1 Wavefunctions 

The first check of the system's dependability was an analysis of the lowest 
order wavefunction of the system, the first subband. For an infinite square 
well effective potential, the analytic solutions go like sin(nx) where n is an 
integer corresponding to the subband number. 

Solving for the wavefunction is not a trivial task, however, as one must 
pick it out from among the 2M randomly sorted, eigenvectors, U+ and U 
in the system. Fortunately, the code had already sorted the eigenvectors into 
conducting and non-conducting modes. Since only the conducting modes' 
eigenvectors represented the real solutions to the Schrodinger equation, only 
the current-carrying modes needed to be sifted. 

This sorting of the current carrying modes was achieved using the re- 
lation e tk = a. This derives from the translational invariance of the Bloch 
Function: 



(4.8) 



where t/j n is wavefunction one lattice point to the left of ipn+i- Since an 
infinite square well potential has a translationally-invariant potential in the 



x-direction, the eigenvalues a that solve the homogeneous equation (3.28), 
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are identified as e tk . 

The largest value of k corresponds to the lowest order wavefunction, 
because the dispersion relation for each subband is a parabola, and the 
subbands are placed one above the other. Thus, for a given intersection 
of E, the Fermi energy, with the dispersion curves, the highest value of k 
will be on the lowest energy subband. Once the lowest order wavefunction 
for both positive and negative momentum is identified, they can be plotted 
across the width of the sample (i.e. in the y direction). 

At zero magnetic field, these two lowest order wavefunctions, simple 
sin 2 (a;) waves (since these are the lowest order eigenvector solutions under 
an infinite square well potential), should be indistinguishable, and indeed 



they are (see Figure 4.1 where the positive momentum wavefunction is scaled 
up by a factor of 1.1 for clarity). As the magnetic field is increased, it is 
expected that the Lorentz force should begin to impact the wavefunctions, 
pushing them to opposite ends of the channel (since they move in opposite 



directions). This splitting is precisely what was found (see Figures 4.2 4.4). 



Higher values of magnetic field meant the wavefunctions were pushed further 



towards the edges of the system (compare Fig. 4.3 with Fig. 4.4). 



4.3.1.2 Real Bands 

The next test conducted to make certain the system, particularly the eigen- 
values and eigenvectors, behaved as expected was to look at the real band 
dispersion relations of the system. Since the value of k (specifically k x ) is 
buried inside the wavefunction, one method to extract its value is to take 
advantage of the Bloch function's result, e lk = a. 

As mentioned above, the dispersion relation between E and k is parabolic 
for a ID system. This is because solving the Schrodinger Equation for this 
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Figure 4.1: Wavefunctions generated with numerical Green's function code 
at zero magnetic field. Red (upper trace) represents the lowest- energy wave- 
function moving to the right, blue (lower trace) the lowest-energy wavefunc- 
tion moving left. The plots run across a slice, from y=-M/2 on the left 
to y=M/2 on the right. The positive wavefunction is scaled up so that the 
wavefunctions do not sit right atop one another. 




y (lattice units) 



Figure 4.2: Wavefunctions at B=0.1 T. The Lorentz force pushes the wave- 
functions moving in opposite directions against opposite walls. 




y (lattice units) 

Figure 4.3: Wavefunctions at B=0.3 T. 
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Figure 4.4: Wavefunctions at B=l T. 
situation one obtains: 

h 2 k 2 

E kx , n = E 0>n + — (4.9) 
2m* 

where £b,n = (n+ \)^^ is the subband energy. 

At B = 0, one therefore expects plots of each subband to look parabolic. 
Furthermore, for an infinite square well potential, it is well known that the 
solutions are spaced apart like n 2 , and hence at B = 0, one expects the 
parabolically shaped subbands to be spaced increasingly far apart. The 
numerical Green's Function model confirmed these expected results. 

For a system with non-zero, perpendicular magnetic field given by the 
Landau gauge, however, the Lorentz force creates an effective confining po- 
tential given by: 



1 eB fik 

^e(y) = 2 m * w c(y - yo) 2 where uj c = — & y = (4.10) 
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As the magnetic field increases, the confinement due to this potential in- 
creases. As B rises, the magnetic confinement potential will have a greater 
impact on the subband energies than the infinite square well potential will. 
This is especially true for the lower bands, as the potential arising from 
the magnetic field is subband independent, while the higher the subband 
number the greater the subband energy due to the square well potential. 
Since the energy spacings resulting from a harmonic oscillator potential 



(4.10) are evenly separated, one would expect to see the spacing between 
subbands change from increasing spacing (like n 2 ) to even spacing as mag- 
netic field rises. Moreover, the spacing between lower order subbands should 
become even first. This result was exactly what was found using the numer- 
ical Green's function model to plot real bands at various magnetic fields. 

An additional feature worth confirming is the shape of these dispersion 
curves. At B = 0, one expects to see strictly parabolic-shaped bands. That 
is, if k is made to be non-zero then the energy must increase parabolically. 
This is not the case for B^O. In this scenario, because the cyclotron orbits 
fit within the square well potential boundaries, even if k is increased, the 
energy does not necessarily increase. That is, the energies become disper- 
sionless once B is sufficiently large. Therefore, the plots of the real bands 
become flat bottomed, not parabolic shaped, at sufficiently large B. The 
program again confirmed these expected results. 



4.3.1.3 Comparison with Analytic Evaluation 

As a third test of the model, local DOS plots using the Green's function code 
were compared to plots generated using the analytic solutions for the local 
DOS. For a one-dimensional quantum system, the system's carrier density 
is given by: 
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Tin = -^•y ,/ 2m*(e - e , n ) for e , n < e (4.11) 

where eo, n is the subband energy. The density of states across a slice, p(m), 
can be given not only by the method derived above, but also by taking the 
derivative of the carrier density of each subband multiplied by the probabil- 
ity map of each wavefunction: 

. , . ,9,dni . ,o,dno , ,o,dn 

rf™> = itfi-*+h*i-5r+---+ <" 2) 

where p is the number of conducting subbands. This equation can be seen 
as equivalent to the general definition of the density of states presented in 



(4.3) : 



p{E) = Y,^-e n )H 2 \ (4.13) 

n 

Given the parameters of the eigenvectors, the DOS, for positive momen- 
tum only, is equivalent to: 

*n)™± , H±faj (4.14) 

where m* should not be confused with the lattice position, m, and U^_[m, I] 
is the probability of finding the Zth conducting positive eigenvector at lattice 
point m on any slice (assuming a translationally-invariant potential). Equa- 
tion ( |4.14| ) is only half the summation for local density of states at point m, 
as U?.[Z,m] must also be taken into account. By calculating p(m) at every 
lattice point on a slice, it can be plotted as a function of m across the slice. 

As both e and B are varied, it is found that the analytical expression 
and the Green's functions provide identical plots. This correspondence be- 
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tween values for the local density of states found using Green's functions 
and the local DOS calculated using an exact solution, evince yet again that 
numerical Green's functions derived from a tight-binding Hamiltonian pro- 
vide an extremely accurate model. The specific behavior of the local density 
of states along a slice is analyzed in the following two sections. 

4.3.2 Behavior of the Density of States at Zero Magnetic 
Field 

So long as the potential is translationally invariant in the x-direction, then 
the density of states plots across a slice should be identical for every slice. 
When surface gates are attached to the sample, however, this translational 
invariance is broken and the density of states will vary for different slices. 
Therefore, without gates, it does not matter which slice is considered. This 



is confirmed both by 2-D plots (see Figure 4.9) and by mapping the DOS 
across a slice for many slices and subtracting them to show they are in fact 
identical. 

Considering first the case of no applied magnetic field, the system is 
examined as the Fermi energy increases. One observes that as E grows, the 
number of bands grows as well. This is to be expected: at a higher Fermi 
energy there ought to be more available eigenenergies. One further sees that 
at the point where E attains a value just greater than an eigenenergy — and 
hence a band appears — the density of states across the sample becomes very 
large. This result is exactly as is to be expected for a 1DEG given that the 



DOS is the derivative of the carrier density (see (4.11)): 



p(E)=2 £ ^p(£-£ 0in )-§ (4.15) 

n,E , n <E 
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Given the form of this equation, there should be sharp peaks in p{E) where 
E — Eon is a minimum, that is, where the subband energy is just less than 



the Fermi energy E. This is precisely what was found (see Figure 4.5) 

The physical reason underlying this peak behavior results from the den- 
sity of states being directly related to the sum of the wavefunctions squared 



(4.12). These DOS plots should exhibit the shape of the highest subband's 
wavefunction squared, since the highest subband will have the highest value 
of Eo tTl and thus the highest value of (E — -Eb,n)~ 5 and therefore dominates. 
Thus if, for example, the third subband's energy, £"0,3 is just less than E, 
there should be three sharp peaks with the troughs very near zero. The 



plots using Green's functions confirm this expectation (Fig. 4.5). Since no 
magnetic field is present in this case, the positive and negative momenta 
wavefunctions overlap. This is not true under the presence of a finite mag- 
netic field. 

In the opposite case, where E is not near a subband's energy, one sees 



that no wavefunction dominates (see Figure 4.6). Therefore, the density 
of states slice plots should show a sum of all the subbands' wavefunctions 
squared. For the case of the infinite square well potential, the density of 
states plot should look like J2 n sm 2 (nx) graphs, since the solutions for a 
system in an infinite square well are sin 2 (nx) solutions. A sum of such 
solutions has troughs that do not come near the x-axis like a single sin 2 (x) 
plot; instead, both the peaks and the troughs should be well elevated from 
the zero-density level. As expected, using numerical Green's functions, the 
shape of a sum-of-sines-squared plot is readily apparent in the DOS plot 



(see Figure 4.7). 

In an infinite square-well potential, the results were identical for every 
slice of the system as can be seen from both a 2-D plot of the density of 
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y (lattice units) 

Figure 4.5: Local DOS at Fermi energy just greater than that of the third 
subband. Generated using numerical Green's Functions. 




y (lattice units) 



Figure 4.6: Local DOS at Fermi energy well above the energy of the third 
subband, but below that of the fourth. Generated using numerical Green's 
Functions. 
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y (lattice units) 

Figure 4.7: Local DOS at Fermi energy well above the energy of the third sub- 
band, but below that of the fourth, compared with plot of sin 2 (x) + sin 2 (2x) + 
sin 2 (3x) ; (green trace). 
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DOS Across Slices 




10 20 30 40 50 60 70 80 90 

x (lattice slices) 

Figure 4.9: Local Density of States in two- dimensions for an energy of three 
subbands at zero field. Black regions represent location of high DOS, yellow 
low DOS. Calculated using numerical Green's functions model. Note the per- 
fect translational symmetry because of the translationally- invariant potential 
used. 



states (Figure 4.9), or from taking 1-D slices from the across the system. 



This latter method was carried out for an iV=50, iV=100, iV=200, and 
7V=500 systems with ten different slices from throughout the sample — in 
some calculations evenly spaced and in others unevenly spaced — and each 
DOS was found to be identical (results are not reproduced as they simply 
show a single DOS since every slice's DOS overlaps perfectly). In short, 
for each of the density of states plots, the theoretically expected densities, 
peaks, and number of subbands across all values of E were generated using 
the numerical Green's functions method. 



4.3.3 Density of States at Finite Magnetic Field 

When a magnetic field is applied, the physics of the system becomes much 
more interesting. What happens to the lowest order wavefunctions upon 
application of a magnetic field has already been examined: the Lorentz 
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force causes the wavefunctions moving in opposite direction to move against 
opposite walls of the system. Now the more general situation of the density 
of states is considered. This involves the behavior of multiple wavefunctions 



and the value of the Fermi energy (4.12). 

At arbitrary Fermi energy, as the magnetic field is increased, the number 
of subbands decreases. This results from the energy solutions to Schrodinger's 



equation under a magnetic field (4.9), in which Eq iTI ~ B. As the increase in 
magnetic field causes the subband energies to grow, then at a given Fermi 
energy, there are fewer subbands. Thus at B=0, the number of subbands 
is at a maximum; this number gradually decreases as B grows until all the 
subbands are depopulated. 

Yet, in examining the local DOS plots, one obtains results that offer 
greater detail than subband number alone. For while the number of sub- 
bands decreases as B grows, the number of peaks in the local DOS plot 
increases. If one depopulates a subband at low enough B-field, there will 
be a drop in the number of peaks in the DOS. Yet, if one were to continue 
to increase the magnetic field, a large new bump would appear in the den- 



sity of states calculation (see Figure 4.10). At first glance, this bump looks 
identical to a new conducting subband entering the system; however, this 
is not the case. For at higher magnetic-fields, beginning most noticeably 
around 0.5 T, the Lorentz force begins to separate all the wavefunctions, 
pushing positive and negative conducting modes to opposite edges of the 
sample. This means that the local density of states plots contain not simply 
a sum of wavefunctions forming a J2 n s'm 2 (nx) graph as before, but a sum 
of positive momenta wavefunctions and a sum of negative momenta wave- 
functions each moving towards opposite edges of the system. As the B-field 
peels apart these wavefunctions, extra bumps, from new overlaps, appear in 
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Figure 4.10: As the magnetic field is increased, one subband depopulates 
at B = 0.3T (lower trace), bringing the total number of subbands to two. 
Increase the magnetic field further and a new hump appears, the result of the 
wavefunctions being pushed apart by the Lorentz force, and the second peak 
is the positive and negative momenta of the n=2 wavefunction overlapping; 
calculated at B = IT (upper trace). 

the DOS plot. 

4.3.4 Imaginary Band Structure 

The local DOS program can also be used to study imaginary band structure. 
Imaginary band structure displays interesting and not entirely understood 
behavior relating to the electronic properties of a quantum system and has 
been investigated, along with the more general complex band structure, in 
a variety of studies JS01EI1E21IM1IM1E3IM1E3EH]- Here, the imaginary 
band structure is derived in much the same way as the real band structure 
was earlier, using the eigenvalues of the transfer matrix. The imaginary 
bands are calculated at the n=0 slice. The difference from the real band 
calculation is that instead of using e tk = a, what is required is: 



e^ ik ) = a — ► k = -lna (4.16) 
The imaginary band structure for a translationally-invariant potential is 
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Figure 4.11: Imaginary band structure, E as a function of Im{k), at B = 0, 
under the presence of no surface gates and assuming an infinite square well 
potential. E on the vertical axis runs from to 10 meV. Calculated using 
numerical Green's functions. 



plotted in Figures 4.11 4.13 At zero magnetic field, these negative paraboli 
are expected given the parabolic dispersion relation between E and k. The 
maxima of these paraboli correspond precisely to the minima of the paraboli 
of the real band structure. As the magnetic field is turned up, however, the 
structure of the subbands becomes more nuanced. A ripple-like effect in the 
system is easily visible at just B= 0.3T. It is also interesting to note the 
asymmetries in the imaginary band structure present at just B=0.7 T that 
become even more pronounced at higher magnetic field. 

Especially intriguing behavior is observed when a surface gate is placed 
close enough to the edge of the system to impact the band structure. A 
similar behavior to the non-gated system for low i?-field is observed, but 
the blending, i.e. the degeneracy, of the bands assumes a different behavior 
in this setup. A tiny rectangular gate of length 6.5 nm and width 26 nm, was 
centered in the middle of a 1-D system of length 65 nm and width 131.3 nm. 
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Figure 4.12: Imaginary band structure at B = 0.3 T, under the presence of 
no surface gates and assuming an infinite square well potential. Calculated 
using Green's functions. 




Figure 4.13: Imaginary band structure at B = 0.7 T, under the presence of 
no surface gates and assuming an infinite square well potential. Calculated 
using Green's functions. 
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Figure 4.14: Imaginary band structure, E as a function of Im{k), at B = 0, 
under the presence of a small central surface gate. E on the vertical axis 
runs from to 10 meV. Calculated using numerical Green's functions. 



The resulting imaginary band behavior can be seen in Figures 4.14 4.16 



Initially, the nature of the points where two imaginary subbands merge 
was investigated, in particular where subbands 3 and 4 or 5 and 6 overlap. 
The local density of states above and below these points of degeneracy was 
plotted, and each 2-D density map was subtracted from the other, but no 
pattern emerged. The density differences between a plot for a Fermi energy 
above the degeneracy and for a plot below the degeneracy could be just as 
prominent as the density differences between two plots both energetically 
above the degeneracy point or both energetically below it. 

What turned out to be the most interesting facet of the imaginary band 
structure was that at high magnetic field the number of conducting modes 
depended upon the bottom portion of the second imaginary subband. Since 
the number of conducting modes depends only upon the real number of sub- 
bands this seemed a strange property. Between the top of the first imaginary 
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Figure 4.15: Imaginary band structure at B = 0.3 T, under the presence of 
a small central surface gate. Calculated using numerical Green's functions. 




link 

Figure 4.16: Imaginary band structure at B = 1 T, under the presence of a 
small central surface gate. Calculated using numerical Green's functions. 
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Figure 4.17: Imaginary band structure at B = 1.2 T, under the presence of a 
small central surface gate. The horizontal lines demarcate the Fermi energies 
as whose values a subband is added (or, in the strange case, subtracted). 
Calculated using numerical Green's functions. 



subband and the bottom of the second subband the system had two conduct- 



ing modes (see Figure 4.17). Above the bottom of the second subband and 
below the top of the second subband, the system had just one conducting 
mode, then returned to two again at the top of second imaginary subband. 
A plot of the real subbands under the influence of the same nearby gate 



revealed the reason for this behavior (see Figure 4.18). The gate distorts the 
first real subband, causing the energy values for small but non-zero k to be 
less than the values of the subband at k=0. As a result, selecting an energy 
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Figure 4.18: Real band structure at B = 1.2 T, under the presence of a 
small central surface gate. The horizontal lines demarcate the Fermi energies 
as whose values a subband is added (or, in the strange case, subtracted). 
Calculated using numerical Green's functions. 



level below the hump at k=0, but above the two lowest minima on either 
side of k=0, yields four intersections, which corresponds to two subbands 
(each with a positive and negative momentum solution). When the energy 
value just above k=0 is reached, there are only two intersections and hence, 
only one subband. This continues until the second subband is reached at 
which point there are again two conducting modes. 

In this way, the imaginary band structure reflects the real band structure. 
The closed-off ellipse of the second imaginary subband reflects the symmetric 
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dips at small k on the real subbands. The reason for these dips in the 
real subband structure is likely the effect of the gate: states with slightly 
positive or negative momenta will be pushed to either side of the gate as 
B is increased. Those states with no momentum (k=0), however, will flow 
directly into the gate even at high B. Thus, the gate causes real subbands to 
appear for small momenta at lower energies than they do for zero momenta, 
resulting in the anomalous gated imaginary band behavior. 

4.4 The Transmission Coefficients Program 

The transmission coefficients calculation program shared a great deal in 
common with the local density of states program; indeed, almost all of the 
initial work is identical. The aim of the program is to calculate the transmis- 
sion, Tij, and reflection, Ru, coefficients, for leads i and j of the Landauer- 
Biittiker formalism. By calculating these values the degree of transport 
through the quantum sample — a crucial characteristic of ID dynamics — is 
determined. 

Like the local DOS program, the transmission program begins with the 
selection of inputs: M, N, a, E, and B. The program is easily adjusted 
to calculate the transmission coefficients as a function of magnetic field as 
is done in Chapter Five. The effective mass used was again that of GaAs, 
0.067e. The constants u, e, and is are used as before, as are V, V^, and Z. 
The Hamiltonian and its dependence upon an external effective potential 
was again implemented. Eigenvalues and eigenvectors of the transfer matrix 
were found and sorted, and Gq ^° was calculated. 

Here the process diverges from that used for calculating the local density 
of states. There, since the local DOS at every lattice point was sought, one 
had to iterate all the way across the system in order to find the DOS for 
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a single slice and therefore there were N total iterations across the system. 
By contrast, finding the transmission coefficients requires iterating across 
the system only once. 

In the local DOS program, the end goal is to find G^2 matrix; in the 
transmission coefficients program, the goal is to find matrices G^ N (of 



(3.23)) and G ^J (of (3.25)), then add on the leads. The latter of these 
matrices is particularly important as it relates the Green's function at one 
end of the system to the other, i.e. transmission. To determine them, G §° 



is inserted into (3.23) as before to find G^\, which, in combination with 
^0 o° yields GqI- The repetition of this process N iterations later gives 



.(AO 

'N,N 



(of (3.23)) and G 



(TV) 

ojv- 



This iterative process is followed by a series of matrix multiplications 
that are carried out in order to string together Green's functions that fully 
describe transmission across the system. The many matrix products are 
then incorporated into an expression that determines the transmission, T{j 
from any one subband i on one side of the system to subband j, on the other. 
The expression is given by the sum of four terms (with Green's functions 
condensed wherever possible): 



1) U+*G+£VG+- 1)JV+1 U+a+VV 

2) U+*a+*G - °°VG+^VG+- jiV+1 VtG+^ +1 U+Vtvt 

3) -U+*a+*G^VG+^VG+^ 1)JV+1 U+a+VtV 

4) -U+*G+^VG+- ljAr+1 VtG+^ +1 U+VVt 

This expression is derived from the work of Appendix B of [53], where 
the lattice form of the transmission coefficients in terms of Green's function 
is derived. An operator related to the current density, K op (n), is introduced: 
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K op {n) = (V ntm \n,m){n + l,m\ - V^ m \n + l,m)(n,m\j (4.18) 

Where V n . m and V^ m are the mth entry of matrices V and V^, respectively, 
calculated at slice n; \n + l,m) is an eigenvector one lattice point to the 
right of the eigenvector at point n, m. K op {n) can be thought of as the 
current that passes between slice n and slice n + 1. Furthermore, K op (n) is 
very similar to the current density operator J p(n) used to link the Green's 
function formalism to the scattering formalism (see [53 for details), and 
therefore can be related to the conductance coefficients, gij, between leads i 
and j. This in turn leads to an expression for the transmission coefficients, 
where Gjj is the Green's function connecting leads i and j across the entire 
system: 



ih 



kj,ab = — z&alKopin^GijKopiri)^ 



where i ^ j 



(4.19) 



Here, ip^ an d V'h are the eigenvectors of subbands a and b moving in the 



positive and negative directions, respectively. Inserting (4.18) into (4.19), 



one obtains four terms, which when expanded give the terms of (4.17). The 



ip a and ipb vectors become the U+ and U_ entries; the 14,m terms correspond 
to V and VT; the eigenstate and Green's function products form the basis 



of the Green's functions multiplications in (4.17). This is the source of the 
transmission program. 



Once (4.17) is applied, then for the case of only two leads, one obtains 



a transmission matrix whose entries are simply Uj. Then, as given by the 



Landauer formula (2.13), one sums over every value of the matrix. This 
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process calculates the total conductance, across all subbands, through the 
system. 

A similar method is applied to find the reflection coefficients. The rel- 
evant expression for calculating the reflection coefficients from lead i back 
into lead i is Ra. It is given by the sum of four terms (with Green's functions 
condensed wherever possible): 



1) U-*G+^ ljAr+1 U+a+VV 

2) U-*a-*G+^VG+- jiV+1 VtG+^ +1 U+Vtvt 

3) -U-*a-*G+^VG+°^ +1 U+a+VtV 

4) -U-*G+- jAr+1 VtG+^ +1 U+VVt 

This expression is derived from [53 , where the reflection coefficients 



equivalent to (4.19) can be found. The result (4.20) can be derived using 
identifications similar to those listed above for the transmission coefficients. 
As a check on the system, the sum of the total transmission and reflection 
coefficients should be equal to the number of total conducting subbands. 
This was found to be true. 

Results from the transmission program for an ordinary system matched 
theoretical expectations. For the case of a translationally-invariant poten- 



tial such as in Fig. 4.9 the transmission coefficients program produces 
three transmitted and zero reflected subbands, as expected. Also as pre- 
dicted, as the magnetic field is increased leading to depopulation, the pro- 
gram calculates that the number of transmitted subbands falls to zero while 
the number of reflected subbands stays fixed at zero. Under the presence 
of a non-translationally invariant potential — for example, a surface gate in 
the middle of the channel — some subbands are reflected. The reflection 
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coefficient will actually decrease at higher magnetic field, both because of 
magnetic deopopulation and because the Lorentz force pushes the subbands 
against the the sample's edges, allowing them to move around the central 
gate and conduct to the opposite side of the sample. 

What is intriguing about the transmission coefficients program, however, 
stems not from its values in various geometries alone, but from comparing 
its results with those from the density of states program. Here, the limits 
of experimental measurement — which can only determine conductance — are 
made evident. The density of states model can provide details of quantum 
systems unseen by direct experiment and previously understood only by 
inference. The next and final chapter considers two of the many future 
directions to which these two programs could be put in combination to 
probe frontiers in low-dimensional quantum physics. 



Chapter 5 



Future Directions 



5.1 Introduction: Integrating the DOS and Trans- 
mission Coefficients Programs 

As the local DOS is a fundamental property of a quantum system and the 
conductance through a ID system is the fundamental limit of studying elec- 
tronic transmission, there are innumerable future investigations in which 
the programs presented in the previous chapter could prove of service. Any 
experiment relating to transport through a narrow constriction (narrow, be- 
cause the wider the system is, the exponentially longer the calculations be- 
come), inclusive of most geometries, attempting to probe conductance, ther- 
mopower, magnetic-field effects, energy eigenstate distributions, and other 
quantum properties could be modeled with these two programs. Below, it 
is examined how the models presented in this thesis could be applied to 
two current realms of inquiry: antidot behavior and zero-dimensional to 
one-dimensional tunneling. 
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Figure 5.1: An Antidot is a potential hump frequently placed in the middle of 
a split-gate ID channel. Circular orbits represent edge states moving around 
the antidot. Image from [Wj. 



5.2 AntiDot Behavior 



An antidot is a bump in the effective potential in a ID channel, producing 



a region from which electrons are excluded (see Figure 5.1). It can be 
formed by placing a very small gate in the middle of a split-gate device. 
Antidots contain many interesting properties and are ripe for the study of 
magnetoconductance, scattering, and tunneling behaviors |69S EOJ El] • They 
have been studied in a variety of situations, often with intriguing results 
concerning Quantum Hall edge states as well as spin properties [721 1751 173]. 

Antidot systems are excellent candidates for analysis by the local DOS 
and transmission coefficients programs and might be especially well deployed 
in a single antidot system described in [76 . Here, an eigenstate moves 
around the antidot. Tunneling is present both across the channel, leading 
to increased transmission, and tranverse to the channel, leading to increased 
reflection. Varying an applied magnetic field causes resonant dips and peaks 
in the transmission and reflection coefficient values. The voltage on the 
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Figure 5.2: Conductance as a function of magnetic field through an antidot 
system. A value of one means one subband (one moving left, one moving 
right) passing through the system. Voltage on antidot=-3.5 mV. Note the 
resonant dips due to increased reflection. 

antidot (and therefore the radius of the antidot) is an essential parameter of 
the system, because changes in the antidot 's size affect the eigenstates and 
therefore the degree of tunneling that is possible. 

Applying the transmission and local DOS programs to calculate the 
transmission coefficients as a function of magnetic field revealed sharp res- 
onant dips at places where tunneling transverse to the channel suddenly 
peaked. These peaks cause transmission across the channel to decrease 



(see Figure 5.2). One would expect making the voltage on the central gate 



more negative (thereby increasing the size of the antidot), would push the 
eigenstates around the dot closer to the ones at the channel's edges. This 
increased proximity would increase vertical tunneling and therefore cause 
an onset of resonant vertical tunneling at lower magnetic field. This is pre- 



cisely what was found (see Figure 5.3 where the resonant dips are apparent 
at lower B-field). 

Results of greater interest emerge when the local density of states pro- 
gram is applied under the same conditions to the same system. This provides 
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Figure 5.3: Conductance as a function of magnetic field through an antidot 
system. Voltage on antidot=-4-5 mV. Note how the more negative voltage 
lowers the magnetic field at which resonant reflection begins to occur because 



there is greater subband overlap at lower B-field than in Figure 5.2 



for a portrait of the sample's transmissive behavior. The decrease in tun- 
neling across the system is evident in the vanishing cross-channel density of 



states at the values of the resonant dips (for example, see Figure 5.4 the 



local DOS at B = .93 T). The edge states at these resonances are not trans- 
mitted across the channel. If one looks at a plot of the local DOS where the 
conductance is nearly one again, however, the local DOS makes plain that 



transmission across the channel is once again present (see Figure 5.5). 

The local DOS was calculated at dozens of magnetic fields in the sys- 
tem. For a iV=200 system of width 250 nm, calculating each local DOS plot 
for a given magnetic field takes under five minutes on a standard desktop 
computer, and the expediency of the code itself could always be improved. 
Applying a whole host of new conditions could be simply done and studies 
of edge states, tunneling, spin properties, and other phenomena could be 
carried out and analyzed rapidly by looking at the transitions of the trans- 
mission coefficients and local DOS across a range of energetic or magnetic 
perturbations. 
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Figure 5.4: Local DOS plot for B = .93 T. Blue represents high local DOS, 
yellow low local DOS. Transmission is imperceptible here. 
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Figure 5.5: Local DOS plot for B = .96 T. Transmission is evident here in 
the blue lines that cross the entire system which are the subbands pushed to 
the edges by the Lorentz force. 
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5.3 0D-1D Tunneling 

A second active area of research where this pair of programs could be de- 
ployed is the investigation and modeling of zero-dimensional to one-dimen- 
sional tunneling systems. Most recently these systems have been probed 
using a surface-acoustic wave (SAW) to create a moving quantum dot — a 
zero-dimensional, fully confined electronic system — that then tunnels into a 
neighboring 1DEG [77]. These tunneling oscillations could prove crucial to 
creating a qubit, the backbone of a quantum computation device. 

To model such a device, one defines a series of surface gates using the 
formalism derived from |57j . In essence, one creates two ID quantum chan- 
nels, side-by-side, with a small, electrostatically-defined potential hump sep- 
arating them; this hump is controlled by a gate, the gate voltage in turn 
controlling the degree of tunneling between channels. In one of the ID chan- 
nels, the potential due to a surface acoustic wave is modeled as a sine wave 
moving in the x-direction and boundaries in the y-direction decaying sharply 
like Fermi-Dirac functions against the edges of the ID channel. In initial 



calculations performed, the potential was approximated as (see Figure 5.6): 



Asln l A/4 -2> . 
rSaw ( e 10{V-Viop) + l)( e -10( y -Jfc o ttom) + 1) 1 J 

where ybottom an d ytop are the electrostatic boundaries of the 1DEG channel 
in which the SAW is present, x m id is the middle slice of the system, A is the 
amplitude of the saw in mV, and A is the wavelength of the SAW, on the 



order of 1 /mi. Using (5.1) as the effective potential, adding in the effect 
of the gates, and then choosing suitable lattice dimensions for the channel, 
the numerical models presented above tap into stores of information about 
the expected theoretical behavior of the resonant tunneling system. For 
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Figure 5.6: Effective Potential for SAW in 0D-1D tunneling system, using 



(5.1). White demarcates low effective potential, light blue medium potential, 



dark blue high effective potential. The SAW is clearly confined to the upper 
channel with the moving quantum dot at the upper channel's center. 

example, one can analyze the distribution of the local density of states, 
identifying potentially anomalous behaviors in wavefunction distribution. 
New gate geometries can be quickly modeled and potential problems in 
transmission immediately identified and calculated. 



The effective potential of (5.1 ) considers the SAW at a single snapshot in 



time, at the moment when it has a potential minimum midway through the 
channel. A time-dependent potential could be incorporated into the Green's 
function model to make it an even more powerful probe of the system. In 
this way, the truly dynamic nature of the quantum dot due to the SAW could 
be modeled. One simple method of probing the local DOS behavior as the 
SAW moves through the system would be to run the calculations for the 
SAW with its potential minimum at several hundred different lattice points, 
creating a "moving" picture of the dynamics of this tunneling system. 
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5.4 Conclusions 

This thesis has presented a model for calculating the local density of states 
and transmission coefficients. The foundation of the program — the itera- 
tion of numerical Green's functions to obtain relevant energetic information 
about the system — has been derived and explicated. The basics of semicon- 
ductor transport in a 1DEG, employing the Landauer-Biittiker formalism, 
have also been elucidated. 

Two numerical programs have been proposed and presented. The ex- 
tremely high degree of accuracy of these programs has been demonstrated 
in several manners. The programs' results, most notably the local DOS pro- 
gram, match analytic results for local density of states, real band structure, 
and waveband nature and behavior with great precision. The local density of 
states program also evinces expected theoretical behavior in terms of num- 
ber of subbands, peak positions, and depopulation under the influence of a 
magnetic field and free of its influence. Furthermore, it has been shown that 
the local DOS program can be used to bring to light interesting difficulties 
in the imaginary band structure and leaves room for others to be investi- 
gated. The local DOS can be used to quickly provide detailed schematics of 
wavefunction behavior under the presence of almost any arbitrary effective 
potential. Together, the transmission program and the local DOS program 
are easily adaptable to the investigation of numerous ongoing experimental 
inquiries, including 0D-1D tunneling and antidot systems. 
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Appendix A 



Derivation of Iterative 



Equations 



The four iterative Green's functions equations for relating two slices of the 
quantum systems are given by: 



Gi n + + ^ +1 = [Z - H n+1 - VtGSvr 1 (A.l) 

Gg +1) = GW + GgVG&^iVtGg (i,j<N) (A.2) 

G&S = GgVGSSSlf! («<JV) (A.3) 

= G^i^VtGg (i<iV) (A.4) 



Equation (A.3) is derived directly from the Dyson equation for the case 



j = n + 1. The first term drops out because the nth Green's function does 

(n) 

not contain an entry for G- n , x . This is because by the nth iteration, there 



is not yet an n + 1 column. Only once G (n+1) is calculated can uA.3\ provide 
the values of the last column of the total Green's function matrix — save the 
very last row's entry. 
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APPENDIX A. DERIVATION OF ITERATIVE EQUATIONS 



84 



This final row, final entry of the Green's function given by (A.l) is de- 



rived from (A. 3). First i,j = n + 1 are plugged into the Dyson equation. 



Second, the result of (A. 3) is substituted in for G^jJ"^ yielding: 



^(n+i) _ {i(n) , p.W vtrWvr( n+1 ' 

^n+l.n+l — ^71+1,71+1 "T ^n+l.n+l v ^"n,n v ^n+l.n+l 



.(«) 



(A.5) 



where is used instead of V because it connects G^"^ n+1 to a lattice slice 



(n) 

to the left, G n ,n- Subtracting the right-most term from (A.5), factoring out 
^n+iji+x on the left-side and multiplying by (G^ lin+1 ) _1 one arrives at: 



,(n+l) 
r n+l,n+l 



((GjJi^-'-VtGWV) 



(A.6) 



which, by the definition of the Green's function reduces to (A.l ). 



Like (A. 3), (A.4) can be derived straightforwardly from the Dyson for- 



mula. If i = n + 1, then as a result of the (j < N) constraint of (A.4), the 



relevant Dyson formulation will be: G^j = G^ , + G^ n+1 VtGjf +1) . 



As with the derivation of (A.3), and for the same reason, the first term 



does not exist. Using (3.22) to substitute for G 



,(n+l) 
'n,j 



and multiplying by 



(G 



n+1, n+1 



) 1 the result is obtained: 



((Git^n 1 - VtGS>V)Gl n+1) 



n+1 j 



The parenthetical term on the left-hand side is (G 



yt G (n) 

(n+1) 



n+l,n+D 



(A.7) 



_1 according 



to (A.6). Inputting this Green's function into (A.7) and multiplying both 



sides by its inverse, returns (A.4). 



The final equation, (A. 2), is derived from (A.4) in a single step. Employ- 



ing (3.22) directly and substituting for the Green's function G^J 1 ] using 



(A.4) gives (A.2). 



